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Abstract 

The  nonlinear  instability  of  a two-dimensional  single  crystal  of  pure  material  grow- 
ing from  an  undercooled  melt  is  studied  both  analytically  and  numerically.  The  qucisi- 
steady  state  approximation  is  used  for  the  thermal  fields  and  the  effects  of  different 
solid  and  liquid  thermal  conductivities  and  isotropic  surface  tension  are  included.  A 
bifurcation  analysis  is  performed  by  calculating  the  instantaneous  value  of  the  funda- 
mental component  of  the  local  normal  growth  speed  for  an  interface  perturbed  by  a 
single  Fourier  shape  component.  The  base  state  is  time  dependent,  and  two  bifurcation 
criteria  are  studied,  the  relative  stabihty  criterion  according  to  which  the  time  deriva- 
tive of  the  ratio  of  the  perturbation  amplitude  to  the  radius  of  the  underlying  circle 
vanishes,  and  the  absolute  stability  criterion  according  to  which  the  time  derivative 
of  the  perturbation  amplitude  vanishes.  Numerically,  the  fundamental  component  of 
the  interfacial  growth  speed  is  found  by  Fourier  analysis  of  the  solution  to  an  integro- 
differential  equation  (obeyed  at  the  interface)  which  gives  the  instantaneous  value  of 
the  local  normal  growth  speed.  Analytically,  a weakly  nonlinear  expansion  technique 
is  used  to  derive  a solvability  condition  at  each  bifurcation.  Our  analytical  and  numer- 
ical results  are  in  very  close  agreement,  and  therefore  mutually  corroborative.  Landau 
coefficients  are  presented  as  a function  of  the  various  dimensionless  parameters  used  in 
the  model.  Almost  all  of  the  bifurcations  are  subcritical. 
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1 Introduction 


In  this  paper,  we  present  a nonlinear  bifurcation  analysis  for  a single  crystal  of  pure 
material  growing  from  an  undercooled  melt  in  two  dimensions.  Our  model  couples  the  effects 
of  isotropic  surface  tension  with  quasi-steady-state  heat  flow,  according  to  which  the  thermal 
fields  satisfy  the  Laplace  equation  throughout  the  bulk  phases.  We  treat  the  case  of  isotropic 
surface  tension  but  allow  the  conductivity  in  the  crystal  to  differ  from  the  conductivity  of  the 
melt.  The  nonlinear  bifurcation  analysis  uses  both  numerical  solutions  to  a fully  nonlinear 
model  of  crystal  growth  as  well  as  analytical  (weakly  nonhnear  expansion)  techniques. 

To  date,  analytical  work  on  nonlinear  stability  theory  in  solidification  has  focused  primar- 
ily on  exploring  cellular  interfaces  that  arise  during  directional  solidification,  and  have  small 
but  finite  amplitudes.  Wollkind  and  Segal[l]  were  among  the  earhest  authors  to  perform  an 
analytical  weakly  nonlinear  analysis  for  directional  solidification  in  two  dimensions.  They  in- 
vestigated the  weakly  nonlinear  interaction  of  a perturbation  of  fixed  wavelength  with  itself. 
Wollkind  and  others[2,  3,  4]  have  continued  this  work,  including  three  dimensional  calcu- 
lations. Wheeler[5],  McFadden  et  al[6]  and  Ungar  et  al[7]  have  also  performed  analytical 
nonlinear  calculations  and  have  derived  self-consistent  evolution  (Landau)  equations  for  the 
perturbation  amplitude  as  part  of  the  stability  analysis.  Wheeler [5]  and  Dee  and  Mathur[8] 
have  treated  the  case  in  which  a band  of  wavelengths  is  considered,  and  McFadden  et  al[9] 
have  included  anisotropic  surface  tension  in  the  weakly  nonlinear  analysis. 

In  this  paper,  we  extend  the  analytical  analysis  used  in  directional  solidification  to  the 
case  of  an  initially  circular  single  crystal  growing  into  an  undercooled  melt.  Furthermore, 
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we  use  a numerical  method  to  compute  the  instantaneous  growth  velocity  of  the  crystal- 
melt  interface  and  compare  these  results  with  those  obtained  by  an  analytical  method.  The 
numerical  (boundary  integral)  method  is  the  same  developed  to  calculate  the  evolving  form 
of  the  crystal-melt  interface  during  the  free  growth  of  a single  crystal[10,  11]. 

In  the  following  sections,  we  first  present  our  model  of  crystal  growth,  including  the 
integral  equation  that  is  used  to  determine  numerically  the  local  normal  growth  speed  along 
the  interface  at  a given  time.  We  then  give  a brief  description  of  the  numerical  method  used 
to  solve  this  integral  equation  and  present  our  numerical  results.  The  analytical  nonhnear 
analysis  follows,  and  a comparison  between  the  analytical  results  and  the  numerical  results 
is  made. 


2 The  Model 

We  consider  a two-dimensional  single  crystal  of  pure  material,  growing  at  the  expense  of  a 
surrounding  liquid  phase  due  to  the  presence  of  an  isothermal  heat  sink,  at  temperature  Too, 
located  at  a circular  outer  boundary  of  radius  In  this  model,  lengths  have  been  scaled 
by  the  nucleation  radius  R*  = TmI I L{Tm  — Too)  referred  to  the  surface  tension  7 (assumed 
isotropic),  where  Tm  is  the  bulk  melting  temperature,  and  times  by  r = [R*Y / a^S ^ where 
ai  is  the  thermal  diffusivity  of  the  liquid.  Here  S = PlCl{Tm  ~ Toq)IL  is  the  dimensionless 
supercooling,  in  which  pl  is  the  density  of  the  liquid,  cl  is  its  specific  heat,  and  L is  the  latent 
heat  per  unit  volume.  The  quantity  5 is  assumed  to  be  small,  so  that  the  quasi-steady  state 
approximation  (the  Laplace  equations  in  the  bulk)  is  valid.  Further  details  may  be  found 
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in  previous  references [10,  11].  We  have  assumed  for  the  sake  of  simplicity  that  the  densities 
of  the  two  phases  are  equal  [pi  = ps),  that  the  specific  heats  in  the  two  phases  are  equal 
(c£,  = cs)  and  that  there  is  no  convection  in  the  liquid.  The  governing  equations  in  the  bulk 
are; 

= 0 (in  the  melt)  (1) 

and 

'^^Us  = 0 (in  the  crystal),  (2) 

where  Us,l  = {'^s,L  — Too)/{Tm  — T^o)  is  the  dimensionless  temperature  in  the  crystal  (S)  or 
in  the  melt  (L).  At  the  far  boundary,  Ul  = 0 and  at  the  crystal-melt  interface: 

Us,L  = 1 - A",  and  (3) 

Vn  = [-VUl  + PVUs]  ■ n,  (4) 

where  K is  the  (dimensionless)  local  curvature,  V)v  is  the  local  normal  growth  speed,  ^ = 
ks/ki  is  the  ratio  of  the  thermal  conductivities,  and  n is  the  normal  vector  to  the  interface 
pointing  into  the  hquid  phase. 

3 Numerical  Method 

3.1  Integral  Equation  Method 

Eqns.  (1)  through  (4)  may  be  reformulated  into  a single  boundary  integral  equation.  We 
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define  the  variable  a as  the  quotient  of  the  arclength  s to  the  perimeter  St  of  the  interface; 
then,  for  a given  point  on  the  interface,  Qq,  the  result  is, 

Uiioio)  = St  j G{a,ac)VN(a)da 

+ St[1-S]J  U,(a)VG{a,  a„)  ■ n(a)da,  (5) 

in  which  G{a,ao)  is  the  Green’s  function  (see  [11]  for  details)  for  a circle  of  radius  i^oo,  and 
Ui  = 1 — K is  the  interfacial  temperature.  The  solution  to  the  integral  equation  (5)  will  give 
VatIq^o)  (tfie  local  normal  growth  speed)  at  any  instant  of  time.  This  equation  is  solved  by 
numerical  methods  described  previously [10,  11]. 

3.2  Numerical  Results 

We  assume  an  interface  shape  given  by  r = R-\-  A cos{kd)  and  choose  values  of  the  radius 
R,  the  amplitude  A,  the  ratio  of  the  solid  to  hquid  thermal  conductivities  (3,  and  the  dimen- 
sionless parameter  R^o-  We  use  the  numerical  method  referenced  in  the  previous  section  to 
solve  the  integral  equation  (5)  for  the  local  normal  growth  speed  along  the  interface  at  a 
given  time.  Fourier  analysis  of  the  solution  enables  us  to  determine  the  fundamental  Fourier 
component,  V^ik,  A),  of  the  growth  speed  from  the  numerical  solution.  We  investigate  the 
manner  in  which  Viv(^,  A)  changes  as  a function  of  the  amplitude  A of  the  perturbation.  As 
the  amplitude  of  the  shape  perturbation  increases,  the  solution  of  Eqn.  (5)  displays  increas- 
ingly nonlinear  behavior,  affecting  V)v(fc,  A),  and  causing  higher  order  harmonics,  V^ink,  A), 
in  the  velocity  Fourier  spectrum  to  grow. 
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Examples  of  the  nonlinear  behavior  of  the  fundamental  component  of  the  velocity  V^(  fc,  A) 
as  a function  of  the  increase  of  the  amplitude  of  a perturbation  having  a four-fold  symmetry, 


[k  = 4),  are  shown  in  Figures  (1)  and  (2).  In  these  calculations  R < Rcr,  in  which  Rcr  is 
the  critical  radius  at  which  the  fundamental  component  of  the  velocity  vanishes.  Rcr  may 
be  found  by  using  linear  perturbation  theory  and  is  a root  of[12] 


\\n{Rcr/Roo) 


1 - 


A. 

M 


+ 


k{l-k^) 

Rcr 

R 


= 0, 


(6) 


withM=  (l-(|21)“)/(l  + (|2^)“). 


In  Figs.  (1)  and  (2),  each  curve  represents  a specific  value  of  the  crystal’s  mean  radius  R.  In 
Figure  (1),  we  have  taken  the  thermal  conductivities  in  both  phases  to  be  equal,  while  for 
Figure  (2),  (3  = 2.  We  define  the  critical  amplitude  Ani  to  be  the  finite  amplitude  for  which 
= 0,  and  in  Fig.  (3),  A^i  is  plotted  as  a function  of  R under  conditions  given  in 
Fig.  (1).  As  the  departure  of  the  average  crystal  size  below  the  critical  radius  increases,  the 
critical  amplitude  becomes  larger. 

Numerical  analysis  of  the  Fourier  components  of  the  growth  speed  shows  that  the  fun- 
damental Fourier  component  of  the  interface  velocity  obeys,  to  a good  approximation,  the 
equation 

= (7) 

in  which  Ci  may  be  calculated  from  hnear  stability  theory  and  changes  from  negative  to 
positive  as  the  radius  changes  from  below  to  above  the  critical  radius,  and  Ani  denotes  the 
numerical  value  of  the  critical  amplitude  for  a given  radius  R.  Furthermore,  as  the  amplitude 
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A increases,  higher  order  harmonics  A)  will  appear  in  the  Fourier  spectrum  of  the 

velocity.  As  the  higher  order  harmonics  begin  to  grow  large,  the  second  harmonic  of  the 
velocity  field  is  a hnear  function  of  while  the  third  harmonic  of  the  velocity  is  a linear 
function  of  A^  and  so  on.  Thus,  VN{nk,  A)  = C„A"  for  A sufficiently  small.  We  show 
an  example  of  these  results  in  Table  (1).  In  the  next  section,  we  use  an  analytical  weakly 
nonlinear  analysis  to  reproduce  the  numerically  predicted  nonlinear  bifurcation,  and  to  gain 
further  insight  into  the  nature  of  the  bifurcation. 

4 Weakly  Nonlinear  Analysis 

The  analytical  nonlinear  stability  analysis  presented  in  this  work  differs  from  previous 
nonlinear  analyses  for  the  following  reasons.  First,  the  geometry  of  the  unperturbed  single 
crystal  (the  shape  about  which  the  nonlinear  expansion  is  performed)  is  circular,  not  planar. 
Second,  the  base  state  is  neither  quiescent,  steady-state  nor  time  periodic.  Our  analysis 
therefore  results  in  instantaneous  stability  conditions  because  during  free  growth,  the  cou- 
pling of  all  growing  Fourier  components  acts  to  change  the  crystal  shape  in  the  next  instant 
in  time. 

In  order  to  investigate  finite  amphtude  effects,  we  assume  that  a circular  crystal  is  per- 
turbed by  a single  Fourier  component  of  integer  wavenumber  k,  so  that  r(9)  = R-\-  Acos{k9). 
We  fix  the  wavenumber  of  the  perturbation  and  investigate  the  weakly  nonlinear  interaction 
of  a perturbation  of  given  wavelength  with  itself.  Our  strategy  is  to  fix  the  radius,  R,  and 
find  the  amplitude  A^i  for  which  the  fundamental  component  of  the  interfacial  velocity  is 
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zero.  We  assume  that  the  mean  radius  of  the  crystal,  the  shape  perturbation  amplitude,  the 
liquid  and  solid  temperatures  and  the  velocity  of  the  interface  may  be  expanded  in  powers 
of  a small  parameter  e. 


R{e)  = 

n=0 


n=0 


r(«,£)  = Y,  + E 'V-i  = Y 


n=:l 


n=0 


(8) 

(9) 


U{r,  d,  e)  = U^^\r)  + 9)  + 0)  + ^)  + • • • , and  (10) 

Vn{0.  e)  = + eVf:}\e)  + + e^V^^\0)  + • • • (11) 

in  which  cos(fc^),  and  poi  = 0.  The  velocity  Vjv  is  composed  of  Fourier 

components  at  each  order  of  e.  We  substitute  Eqns.  (8)  to  (11)  into  Eqns.  (1)  to  (4),  expand 
all  boundary  conditions  in  a series  about  the  unperturbed  circular  radius,  R^^\  for  which 


e = 0,  and  obtain  a system  of  equations  at  each  order  of  the  expansion  parameter  e.  For  the 
zero  order  problem  we  have 


dr'2  Vr/  dr 

dr'2 


= 0,  for  r < R^°\  with  Us^\r  = 0)  < oo,  and 


(1  \ dU^ 

-J  ^ = 0,  for  R^^^  < r < Roo,  with  = R^)  = 0, 


in  the  bulk  phases,  and 


C/i"’  = 1 - 


R(0)' 


!7[°>  = 1 - and 

^ iJ(o) 


(12) 

(13) 

(14) 
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Vi“>  = + /JDCrf ) . 


at  r = D denotes  the  derivative  of  a function  with  respect  to  r.  For  all  other  orders, 

e",  the  equations  in  the  bulk  are 


= 0,  for  r < with  Us^\r  = 0)  < oo,  and 

= 0,  for  < r < i^oo,  with  Ui"\r  = R^)  = 0, 
and  the  boundary  conditions  are 

c/r + z*"'  (oc/f)  - = c„, 

[/[">  + Z'")  (ZJCf)  - = D„,  and 

- D (-ci"'  + + Z*"*  (D2(C/i‘”  - I3U^S^))  = E„, 


(15) 


(16) 


(17) 


in  which  all  functions  are  evaluated  at  the  unperturbed  radius  R^^K  The  subscripts  6 refer 
to  derivatives  with  respect  to  the  polar  angle  9.  At  any  order,  say  0(e”),  the  inhomogeneous 
terms  in  Eqns.  (17),  C„,  Dn  and  depend  only  on  the  solutions  at  0(e”“^)  or  lower. 
These  inhomogeneous  terms  are  listed  in  the  Appendix.  The  set  of  equations  at  each  order 
in  e is  linear,  and  our  solution  procedure  is  to  solve  this  set  of  equations  up  to  order  e^. 


4.1  Zero  Order  (Unperturbed)  Solution 

The  zero  order  (e°)  problem  gives  the  unperturbed  circular  solution.  The  growth  velocity  is 
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and  the  solutions  for  the  unperturbed  thermal  fields  are 


= 1 - and  (19) 

= -Bo  ln(-^), 

Jtoo 

in  which 


4.2  First  Order  Solution 

In  order  to  find  a solution  to  the  set  of  differential  equations  at  first  order  in  e,  we  assume 
solutions  of  the  form 

+ Bi;iV)cos(M),  (21) 

+ p\\  cos{k9)  and 
vj^'>  = VP  + cos(k»). 

Analysis  of  the  boundary  conditions  (Eqns.  (17))  shows  that  the  differential  equations  and 
the  boundary  conditions  are  satisfied  identically  if  we  set  = 0,  in  which  case 
and  also  vanish.  In  order  to  find  the  critical  amphtude,  we  must  require  the  first  order 
contribution  to  the  fundamental  component  of  the  local  normal  growth  speed  to  vanish. 
This  condition  reproduces  the  marginal  stabihty  results  of  hnear  theory.  Setting 
Ui^\  and  to  zero,  Eqns.  (21)  are  substituted  into  Eqns.  (15)  and  (16),  from 

which  the  solutions 
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and 


(22) 


are  determined.  The  unknown  solution  constants,  and  pn,  satisfy  the  matrix 

equation, 


Pn^ii  — 0, 


(23) 


\ 


.(24) 


in  which  .Yu  = (^5^^  pn)^,  and  Pn  is  the  matrix 

(^(0))*=  0 - ((1 -P)/(P(o))2) 

0 1/(P'')  (1  - (P(o)/Poc)'*'')  (Po/P(°^  - (1  - P)/(P(°0^) 

-^A:(P(o))'=-i  -A:/(P^+i)  (l  + - [Bo/{R^^^)^) 

In  order  to  obtain  non-trivial  solutions  to  the  homogeneous  system,  the  determinant  of  Pn 
must  vanish,  thereby  determining  the  marginal  stability  condition  and  defining  the  corre- 
sponding marginal  stability  radius  = Per-  The  marginal  stability  condition  has  already 
been  given  in  Eqn.  (6).  The  solutions  to  the  first  order  problem  may  be  written  in  terms  of 
a single  undetermined  coefficient  pn  in  the  form: 


^ ^ cos(fc^) 


(25) 


in  the  crystal,  and 


-1 


(26) 


Rt 


R^ 

•^^00 


Pii  cos(A:^) 


in  the  melt. 
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4.3  Second  Order  Solution 


Expanding  the  right  hand  sides  of  the  boundary  conditions,  Eqns.  (17),  gives 

C2  = C20  + C22  cos(2A;^),  (27) 

D2  ~ ■7^20  “I”  7?22  cos(2fc^)  and 
E2  = E20  "I"  E22  cos{2k6) 

in  which  C20,  C22,  -D20,  D22,  E20  and  E22  are  given  in  the  Appendix.  Since  the  second  order 
equations  are  linear,  we  may  assume  solutions  of  the  form 

= UiT(r)  + Cg'(r)cos(2M),  (28) 

and 

^ ym  ^ ^^22)  cos(2^^). 

Equations  (28)  are  substituted  into  Laplace’s  equation,  and  each  Fourier  component  of  the 
solution  of  the  system  of  linear  equations  is  solved  independently.  The  solutions  are 

CrV)  = (29) 

Cr('-)  = -Brin(£-) 

and 

CfV)  = (30) 
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(r) 


Ri^)' 


Eqns.  (29)  and  (30)  are  then  substituted  into  the  interfacial  conditions  whereupon 
and  are  found  from  the  solution  of 


a(^o) 

* ■■rT 


— ^*205 


In 


Rr 


Rr 


and 


d(20) 

y(2o)  + _ R(^) 


Rr 


A 

Rlr 


— ■E'20- 


(31) 


At  this  stage,  we  appear  to  have  six  unknowns,  R^^\  Pii  and  e,  but  only 

three  equations.  A degree  of  freedom,  however,  has  been  caused  by  the  introduction  of  the 
expansion  parameter  e,  and  this  allows  for  an  arbitrary  but  convenient  choice  for  e.  Thus,  to 
second  order  in  e,  only  the  quantities  e^R^^^  = R — Rcr  and  epu  = A have  physical  meaning. 
For  example,  if  we  choose  e = {Rcr  — R)^^^  (or  equivalently  R^"^^  = —1)  in  which  R < Rcr 
is  the  actual  radius  of  the  crystal  (which  would  be  an  appropriate  choice  for  a subcritical 
instability)  the  nonlinear  analysis  could  then  be  continued  to  third  order  to  calculate  the 
value  of  the  critical  amplitude  A = A„/  as  a function  of  i?  = Rni  at  the  nonlinear  bifurcation 
point.  An  alternative  strategy  is  to  choose  e = A (or  equivalently  pu  = 1)  and  then  to 
calculate  the  value  oi  R — Rcr  = Rni  — Rcr  as  a function  of  A„/  at  the  bifurcation  point 
from  the  nonlinear  analysis.  The  sign  of  Rni  — Rcr  indicates  whether  the  bifurcation  is 
“supercritical”  or  “subcritical”.  For  the  moment  we  choose  to  retain  flexibility  and  express 
the  solutions  at  second  order  as  functions  of  the  two  expansion  parameters,  R^^^  and  pn. 
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With  the  aid  of  Eqns.  (31),  the  solutions  and  and  the  velocity  component 

may  be  written  in  the  form 


[/'“>  = Coo  + 


7^(2) 


(32) 


= B20- 


A 

Rrr 


\n{r/Roo) 


i?2 


and 


y^20)  ^ ^ ^(2) 


J?1 


ln(7i’cr/7?oo) 

(|^“  :^))  i^cr-HRcr/Roc)) 


-1 


The  solution  constants  A^g^\  and  are  determined  by  using  the  boundary 

conditions  corresponding  to  the  cos(2A:^)  problem  and  lead  to  the  matrix  equation 


P22X22  — R 


22 


(33) 


in  which  X22  = {A^s^\  ^ R22  = (^*22, -D22, -£^22)^  and  P22  is  the  matrix 

^ (P0)2^  0 0 ^ 

0 1/(P'^)2  (1  - (p(0)/p^)4fc)  0 • 

By  solving  Eqn.  (33),  we  determine  the  solutions 


and 


(34) 


(35) 
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We  emphasize  that  although  the  interface  shape  consists  of  only  a single  Fourier  component 
(proportional  to  cos(^’^)),  there  exists  a Fourier  component  of  the  velocity  proportional  to 
cos(2^’^).  We  may  not  assume  that  this  component  of  the  velocity  is  zero,  for  then  the  analysis 
leads  to  an  overdetermined  system  of  equations.  Therefore,  the  nonlinear  analysis  gives  only 
an  instantaneous  condition  for  the  disappearance  of  the  fundamental  Fourier  component  of 
the  velocity.  Since  this  second  order  solution  does  not  determine  a relationship  between  p\\ 
and  (or  equivalently  between  Ani  and  Rni)  we  proceed  to  the  third  order  problem. 

4.4  Third  Order  Solution 

We  expand  the  right  hand  sides  of  the  inhomogeneous  third  order  boundary  conditions  (17) 
to  find 

C'a  = C31  cos(fc^)  + C33  cos(3A:^),  (36) 

£>3  = £>31  cos{kB)  + D33  cos(3A:^), 

£*3  = £31  cos{k9)  + E33Cos{Zkd), 

From  the  Appendix  we  see  that  C31,  £31  and  £31  are  of  the  form 

C31  = C313P11  + Csiipii  (37) 

D31  = £313/9^1  + D3UP11 
E31  = £313/511  + £311/^11 
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(The  remaining  coefficients  C33,  £>33  and  are  not  needed  for  the  analysis.)  Consistent 
with  Eqns.  (36),  we  assume  solutions  at  third  order  of  the  form 

UfX(r,e)  = u!;%>(r)cos(M)  + [/^“>(r)cos(3M),  (38) 

= P3i  cos(kO),  and 

1/^3)  ^ y^31)  ^ y^33)  cos(3A:^). 

These  are  substituted  into  the  Laplace  equations  and  the  resulting  solutions  are  then  applied 
to  the  third  order  boundary  conditions  to  give 

PiiXzi  = Rzu  (39) 

for  the  cos{k9)  component,  for  which  we  seek  the  unknown  vector  X31  = pzi)^ ■ 

In  Eqn.  (39),  the  right  hand  side  R^i  = ((^31,  £>31,  E3i)^;  furthermore  has  been  set 

to  zero  because  we  are  seeking  the  marginal  stability  condition.  The  determinant  of  Pn 
vanishes  (because  the  radius  has  been  set  at  Per  in  the  first  order  problem);  thus,  in  order 
for  a solution  X31  to  exist,  the  solvability  condition 

{R3uy)=0  (40) 

must  be  obeyed  (in  which  the  parentheses  denote  the  inner  product  of  P31  and  y)  for  all 
vectors  y satisfying  Pfjy  = 0 where  P-^  is  the  adjoint  of  the  matrix  Pn.  Evaluation  of  Eqn. 
(40)  leads  to  the  solvability  condition 

+ (if)  + (^) 
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Substitution  of  Eqn.  (37)  for  C31,  D31  and  E31  into  Eqn.  (41)  leads  to  the  following  system 
of  equations  relating  pn  and  e: 

np?,  - = 0,  (42) 

■^nl  = ^Plli 

R„i  = R,,  + e^R(^K 

The  coefficients  T and  fl  are  given  in  the  Appendix.  Eqns.  (42)  give, 

R„,  - R„  = Al,^  (43) 

which  is  the  central  analytical  result.  In  the  following  section  we  compare  Eqn.  (43)  with 
our  numerical  results  and  characterize  the  nature  of  the  bifurcation  over  a range  of  material 
parameters. 

5 Analytical  Results 

5.1  Comparison  with  Numerical  Results 

To  verify  the  accuracy  of  our  calculations,  we  compare  the  values  of  the  subcritical 
amplitude  Ani  for  given  R^i  determined  by  the  numerical  solution  of  Eqn.  (5)  with  the  value 
found  by  using  the  weakly  nonlinear  stability  analysis.  In  Fig.  (4),  we  examine  the  case  for 
which  the  ratio  of  the  thermal  conductivities  is  unity  and  k = 4.  We  plot  ln(  A„/)  as  a function 
of  1/2  ln{Rcr  — Rni)  and  observe  that  the  numerical  calculations  (denoted  by  the  open  boxes) 
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lie  directly  on  top  of  the  straight  line  of  unit  slope  representing  the  analytical  calculation. 
Our  numerical  and  analytical  results  are  therefore,  mutually  corroborative.  In  Table  (2),  a 
further  comparison  between  the  values  of  the  numerical  and  the  analytical  calculations  is 
presented.  In  Tables  (3)  and  (4),  we  hst  the  values  of  Ani/{Rcr  — RniY^^  corresponding  to  the 
point  of  subcritical  instability  as  the  value  of  the  dimensionless  parameter  R^  increases  for 
fixed  /?,  and  as  the  value  of  the  ratio  of  the  thermal  conductivities  /?  increases  for  fixed  i?oo- 
Table  (3)  indicates  a very  weak  (essentially  logarithmic)  dependence  on  Roo  as  expected; 
Table  (4)  indicates  that  as  (3  decreases,  the  bifurcation  becomes  more  strongly  subcritical. 


5.2  Relative  and  Absolute  Stability  Criteria 


We  have  computed  the  stability  condition  for  which  the  fundamental  component  of  the 
growth  speed  vanishes,  sometimes  called  absolute  stability [12].  (This  is  not  to  be  confused 
with  the  Mullins  and  Sekerka  usage  of  absolute  stability [13].  ) Insofar  as  shape  changes  are 
concerned,  however,  this  is  not  the  most  meaningful  criterion;  even  though  a perturbation 
may  be  growing,  A/R  may  be  decaying.  A more  meaningful  criterion  is  that  of  relative 
stability,  the  point  at  which  the  time  derivative  of  the  ratio  of  the  perturbation  amplitude 
to  the  underlying  crystal  radius  vanishes.  To  third  order  in  e this  condition  is  given  by, 

— f ^ = i ^ / {epu  + eV3i)(Vff^  + 

dt\R)  y i?(o)  + e2i2(2)  j (i?(o) +e2i?(2))2 

Eqn.  (44)  requires  that  at  the  point  of  relative  instability. 


72(0) 


and 


(45) 
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(3.)  _ (VrVn  + 


— 


N 


Ri^)  (i?(0))2 

The  algebraic  differences  between  the  relative  stabihty  point  and  the  absolute  stability 
point  are  the  following;  at  first  order,  the  homogeneous  boundary  conditions  may  be  written 


PnXn  = 0, 

where  is  the  matrix 

/ 


(46) 


\ 


v 


0 0 0 

Ai  + 0 0 0 • (47) 

0 0 ^ 

The  critical  radius  denoting  the  onset  of  relative  instabiHty,  is  now  given  by  the  roots 

of[12] 

fi-{i/R:r) 


\^-h] 

/ 

M 

V Per  ) 

Lm  n 

\MRcr/Poo), 

At  third  order  the  boundary  conditions  are 

P,\X31  = ^31, 

in  which 


= 0. 


(48) 


(49) 


■^31  — + 


v 


(50) 


- vPp,,ir: 

Finally,  the  solvability  condition  for  the  relative  bifurcation  point  is 


/ 


(51) 
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in  which  £^3^  = £31  + pw / puj R*^.  From  now  on  the  superscript  (*) 

notation  will  be  dropped  and  the  stability  criterion  being  considered  will  be  exphcitly  stated. 
We  have  compared  numerical  and  analytical  calculations  of  relative  stabihty  points,  and 
again  we  have  observed  very  close  agreement  between  the  results,  verifying  the  consistency 
of  our  calculations.  Fig.  (5)  gives  a comparison  between  the  values  of  the  relative  and 
the  absolute  critical  radii  given  by  Eqns.  (48)  and  (6),  respectively.  The  relative  stability 
critical  radii  are  larger  than  the  absolute  critical  radii  at  any  wavenumber,  for  fixed  values 
of  /3  and  R^o-  (For  further  details  on  the  condition  of  relative  marginal  stabihty  see  Coriell 
and  Hardy[12].)  Fig.  (5)  is  representative  of  the  difference  between  the  relative  and  absolute 
marginal  stability  radii  for  all  of  the  parameter  values  for  which  we  have  examined  nonlinear 
stability. 

We  investigate  the  subcritical  and  supercritical  nature  of  the  nonlinear  bifurcation  by 
evaluating  the  parameter  (£„/  — £cr)/^^/  as  a function  of  the  conductivity  ratio,  the  pa- 
rameter Roo  and  the  wavenumber  of  the  perturbation.  In  Fig.  (6),  values  of  the  absolute 
bifurcation  parameter  {Rni  — Rcr)l-^\i  are  plotted  as  a function  of  the  wavenumber  for  differ- 
ent values  of  (3  at  fixed  R^o-  The  results  show  that  this  parameter  is  almost  always  negative, 
the  exceptions  being  at  the  largest  values  of  £00  and  /?,  for  wavenumbers  of  A:  = 2 and  k = 3. 
Changing  the  value  of  R^  does  not  change  appreciably  the  features  shown  in  Fig.  (6).  In 
Fig.  (7)  we  show  an  enlarged  view  of  a section  of  Fig.  (6)  that  shows  more  clearly  the  change 
in  sign  of  (£„;  — £cr)/^n/  perturbations  with  wavenumbers  2 and  3.  In  Fig.  (8)  we  plot 
the  parameter  (£„/  — Rcr)IA\i  determined  using  the  relative  stability  criterion.  The  results 
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show  that  according  to  the  relative  stability  criterion,  the  bifurcations  are  nearly  always 
subcritical  except  at  the  wavenumber  k = 3.  The  qualitative  behavior  of  the  bifurcation 
is  very  similar  for  all  cases  considered  for  either  the  relative  theory  or  the  absolute  theory. 
As  the  ratio  of  thermal  conductivities  decreases,  or  as  the  value  of  Roo  decreases,  or  as  the 
wavenumber  increases,  the  bifurcations  become  more  strongly  subcritical. 

6 Conclusions 

In  summary,  we  have  used  numerical  and  analytical  methods  to  calculate  consistently  non- 
linear bifurcations  for  the  fundamental  Fourier  component  of  the  local  normal  growth  speed 
of  a crystal- melt  interface  during  the  free  growth  of  a pure  single  crystal.  The  analytical 
results  correspond  to  an  instantaneous  condition  for  the  nonlinear  stability  of  a given  Fourier 
shape  component,  and  represent  the  first  appHcation  of  weakly  nonlinear  stability  analysis 
to  problems  in  solidification  theory  for  which  the  unperturbed  state  is  not  steady.  Our  nu- 
merical method  is  based  on  a boundary  integral  technique.  The  analytical  analysis  provides 
an  important  check  on  the  numerical  method,  which  can  then  be  used  with  confidence  to 
compute  finite  amplitude  solutions.  Using  the  numerical  solution  technique,  we  have  shown 
that  at  a fixed  radius,  the  fundamental  component,  Vj\f{k,A),  of  the  local  normal  growth 
speed  obeys  a Landau-type  equation.  Furthermore,  the  higher  order  components  of  the 
velocity  are  proportional  to  powers  in  A,  for  small  amplitude  perturbations.  We  find  that 
the  fundamental  Fourier  component  of  the  local  normal  growth  speed  vanishes  at  a critical 
amplitude,  in  agreement  with  results  obtained  from  an  expansion  technique  to  calculate  an- 
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alytically  this  critical  amplitude.  Two  stability  criteria  are  considered;  these  correspond  to 
the  relative  and  the  absolute  stability  criteria.  For  the  absolute  criterion,  our  results  show 
that  the  nonhnear  bifurcations  are  almost  always  subcritical,  except  at  high  conductivity 
ratios  and  high  values  of  Roc  for  wavenumbers  k = 2 and  k = 3.  For  the  relative  stability 
criterion,  the  bifurcations  are  nearly  always  subcritical  except  for  high  values  of  i?oo  and 
wavenumber  k = 3.  Furthermore,  as  the  ratio  of  the  solid  to  hquid  thermal  conductivity 
increases  and  the  outer  boundary  radius  increases,  the  bifurcations  become  less  strongly 
subcritical.  We  also  show  that  for  fixed  thermal  conductivity  ratio  and  Roo,  the  bifurcation 
becomes  more  strongly  subcritical  as  the  wavenumber  increases.  Although  our  results  are 
only  strictly  applicable  to  the  specific  model  that  we  have  treated,  several  broader  inferences 
can  be  made.  First,  the  fact  that  the  results  from  numerical  computations  and  from  weakly 
nonlinear  stability  analyses  are  in  such  good  agreement  is  not  only  comforting  but  suggests 
that  weakly  nonlinear  analyses,  in  general,  make  sense  and  capture  the  most  important  non- 
linearities  near  the  onset  of  instabihty.  This  has  been  appreciated  previously  for  steady  state 
base  states,  but  our  analysis  extends  this  influence  to  non  steady  state  base  states  as  well. 
Second,  almost  all  of  the  bifurcations  that  we  have  examined  are  subcritical  and,  although 
we  cannot  yet  explain  why,  one  is  led  to  suppose  that  this  is  due  to  some  underlying  physical 
mechanism  that  we  should  seek  to  uncover.  One  wonders,  moreover,  if  the  same  results 
(subcritical  bifurcations)  would  be  found  for  a growing  sphere,  and  this  is  currently  being 
investigated.  Finally,  we  note  that  the  nonlinear  analysis  may  be  generalized  to  include  the 
effects  of  linear  isotropic  interfacial  attachment  kinetics;  however,  it  is  not  anticipated  that 
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the  kinetic  effects  will  significantly  alter  the  results  presented  in  this  analysis,  other  than  to 
delay  the  onset  of  morphological  instability  of  the  crystal  to  higher  values  of  R. 
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Appendix 


Definitions: 
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The  coefficients  in  the  expansion  are: 
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Table  1:  Vfs/{nk,A)  as  a function  of  A.  The  numerical  results  show  that  in 
the  regime  A <<  1,  Vi^[nk,  A)  = CnA^ . The  example  shown  in  this  table  is 
for  the  case  given  by  the  curve  labelled  as  (3)  in  Fig.  (1). 


Fundamental  Component 


n 

Cn 

1 

" Cl  = -.452  X 10-« 

2 

C2  = +.363  X 10"^ 

3 

C3  = -.853  X 10"® 

Table  2:  Numerically  and  analytically  determined  values  of  Ani  for  /3  = 1, 
i?oo  = 10®,  Rcr  = 322.56787  and  k = 4. 


Critical  Amplitude 


R 

Ani  (num) 

Ani  (ana) 

322.564 

.730  to  .731 

.730504 

322.563 

.819  to  .820 

.819279 

322.562 

.899  to  .900 

.899333 

Table  3:  Analytically  calculated  values  of  AniURcr  — R-niY^^  as  a function  of 
Roo,  using  the  absolute  stability  criterion  for  ^ = 1. 


Analytical  Analysis 


Rqo 

^nll{Rcr  — RnlY^^ 

10^ 

14.40900 

10® 

16.66126 

5 X 10® 

18.10505 

10^ 

18.69906 

5 X 10^ 

20.02192 

Table  4:  Analytically  calculated  values  of  AniKRcr  — RniY^^  as  a function  of 
(3^  using  the  absolute  stability  criterion  for  Rco  = 10®. 


Analytical  Analysis 


/3 

Ar,l/{Rcr  - RnlY^'^ 

.5 

13.16962 

1.0 

16.66126 

2.0 

22.81829 

Figure  Captions 


Figure  1.  The  fundamental  Fourier  component  of  the  velocity  V0v(^5  as  a function  of  A. 
Here  /?  = 1,  the  value  of  Roo  is  10®,  Rcr  = 322.5678  and  k = 4.  Curves  (1),  (2)  and  (3) 
are  for  radii  R = 322.55,  322.4  and  322.2,  respectively. 

Figure  2.  The  fundamental  Fourier  component  of  the  velocity  VN{k,  A)  as  a function  of  A. 
Here  /?  = 2,  the  value  of  R^o  is  10®,  R^  = 461.8201  and  k = 4.  Curves  (1),  (2)  and  (3) 
are  for  radii  R = 461.6,  461.5  and  461.4,  respectively. 

Figure  3.  Ani  as  a function  of  R plotted  for  the  conditions  given  in  Fig.  (1). 

Figure  4.  A plot  of  ln(A„/)  vs.  l/21n(i7cr  — Rni)-  The  solid  line  is  the  analytical  solution 
given  by  the  weakly  nonlinear  analysis  and  the  open  boxes  are  the  values  from  the 
numerical  calculations.  In  this  case,  /?  = 1 and  R^o  = 10®.  The  slope  is  unity. 

Figure  5.  A comparison  between  the  values  of  the  critical  radii  using  the  absolute  stability 
criterion  and  the  relative  stabihty  criterion.  The  anomalous  case  k = 2 for  the  relative 
stability  criterion  is  omitted  because,  formally,  the  critical  radius  becomes  very  large. 

Figure  6.  A plot  of  {Rni  — Rcr)l-^\i  using  the  absolute  stability  criterion.  The  curves  are 
for  values  of  /5  equal  to  1/2,  1,  2 and  4.  The  value  of  Roo  = 10®. 

Figure  7.  A plot  of  {Rni  — Rcr)/A.'li  using  the  absolute  stability  criterion.  The  curves  are 
for  values  of  /?  equal  to  1,  2 and  4.  The  plots  are  taken  from  a section  of  the  plot  in  Fig. 
(6)  and  reveal  the  positive  values  of  the  parameter  {Rni  — Rcr)/^ni  perturbations 


of  wavenumber  k = 2 and  k = d. 

Figure  8.  A plot  of  {Rni  — Rct)I -A^ni  using  the  relative  stability  criterion.  The  curves  corre- 
spond to  values  of  (3  equal  to  1/2,  1,  2 and  4.  The  value  of  i2oo  = 10*- 
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